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Smearing the gauge links of dynamical configurations removes small scale unphysical vacuum fluctuations and 
thus improves the chiral properties of lattice fermions. Recently we proposed the hypercubic smearing (HYP) 
that improves the flavor symmetry of staggered fermions by an order of magnitude with only minimal distortions 
at small distances. We describe a new algorithm to simulate dynamical HYP fermions based on the standard 
pure gauge overrelaxation and heatbath updates. The algorithm has been used to simulate four and two flavors 
of staggered fermions. Unlike standard dynamical simulation techniques, this algorithm does not loose efficiency 
at small quark masses. 



1. Introduction 

In this article we study an algorithm for simu- 
lating a fermionic system described by the action 

S^Sg{U)-trln[QHV)Q{V)] , (1) 

where Sg{U) is the pure gauge action and Q{V) 
is the fermionic matrix. The gauge connections 
between fermions are smeared links V which are 
constructed deterministically from the dynamical 
thin links U. Since each smeared link is a lo- 
cal combination of a finite number of thin links, 
the system where the fermions couple via smeared 
links is in the same universality class as the sys- 
tem with thin links. 

It has been demonstrated that smeared links im- 
prove flavor symmetry with staggered fermions, 
chiral symmetry with Wilson-type clover fermions 
and they are useful in the construction of an 
overlap Dirac operator. The problem is that 
the standard simulation algorithms for dynami- 
cal fermions (Hybrid Monte Carlo (HMC) or R 
algorithms) involve the computation of the gauge 
force dQ{V)/dU, which is either very compli- 
cated or likely impossible if the smeared links are 
projected onto SU(3). 

2. New algorithm 

We propose for the simulation of the system 
described by eq. (|^) a two step updating 

1) updating of a set of thin links {[/} {U'} 
which satisfies detailed balance for Sg{U) 
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2) accept/reject step with acceptance proba- 
bility 

p . det [QHV')Q{V')] \ 

where V' are the smeared links constructed 
from the updated thin links 

As concerns step 1) our choice for Sg is the Wilson 
plaquette action and the updating is performed 
either with microcanonical overrelaxation or with 
Cabibbo-Marinari heatbath. With overrelaxation 
we update all links within some finite block of the 
lattice (the sequence of the updates in the block 
has to be symmetrized), with heatbath the links 
to be updated are chosen randomly. 
In step 2), we use a stochastic estimator to eval- 
uate the ratio of fermionic determinants 

K;,, = min{l,expA5} (3) 
AS - e[QHV')QiV')~QHv)QiVM, (4) 

where the vector ^ is generated according to the 
probability distribution 

P(Ocxexp{-etQt(r)g(U')e}- (5) 
3. The HYP action 

In the following we describe the application of 
this algorithm to the cases of =4 and nf = 2 
flavors of staggered fermions coupled via hypercu- 
bic (HYP) smeared links. The definition of HYP 
links and their properties are discussed in 
In the natural formulation of four flavors of stag- 
gered quarks, the continuum SU(4) flavor sym- 
metry is broken to U(l) at finite lattice spacing. 
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The pion spectrum has only one true Goldstone 
pion TTG, the other 14 pions remain massive when 
the quark mass goes to zero. The mass sphtting 
between a non- Goldstone pion tt and ttq can be 
parametrized by the quantity 0| 



S2 = 



ml - ml^ 



(6) 



evaluated at rriT^^/mp — 0.55. With a flavor sym- 
metric action S2 — for all the pions tt. We com- 
puted S2 for the two lightest non-Goldstone pi- 
ons (for the notation see on a set of quenched 
8'^ X 24 lattices generated with Wilson pure gauge 
action at /? = 5.7 (a = 0.17fm), using three dif- 
ferent valence quark actions: standard staggered 
action with thin and HYP smeared links and im- 
roved Asqtad action of the MILC collaboration 
The results are shown in table 0. Flavor sym- 



Table 1 

Flavor symmetry breaking for different valence 
quark actions at lattice spacing a — O.lTfm. 



Action/(52 


Thin 


Asqtad 


HYP 




0.594(25) 


0.191(22) 


0.086(14) 




0.72(6) 


0.32(4) 


0.150(24) 



metry violations with the HYP action are reduced 
by an order of magnitude with respect to the thin 
link action and by a factor of two with respect to 
the Asqtad action. 

4. Uf — A HYP staggered fermions 

We consider the matrix Q in eq. ([^) given by 
Q,^j = M,,, =2(am)5,j+Aj(^), (7) 

which describes four flavors of staggered fermions 
coupled via HYP smeared links V . In the ac- 
cept/reject step we rewrite the fermion matrix as 

Q{V) = Qr{V)A{V), (8) 
A{V) = ei^-p[aAD^{V)+a2D^{V)]. (9) 

In this way we achieve that an effective action 

Scff = -2a4Re trD'' - 2a2Rc trD^ (10) 



is removed from the fermion determinant. The 
acceptance probability in eq. (§) becomes Pace = 

min {1, exp[5cff (V^) - S,a{V')] exp A^,} , (11) 

where the stochastic part /S.Sr is computed like in 
eqs. (^^) with Qr instead of Q. The real param- 
eters a4 and 02 are optimized to maximize Pace 
and we use = -0.006, as = -0.18 g. The 
point we would like to emphasize is that this al- 
gorithm is not effective with thin links. Smeared 
links constrain the fluctuations of the stochastic 
estimator exp AS'^ and make the algorithm effl- 
cient. 

In Monte Carlo simulations we separate the mea- 
surements of the observables by A^oR overre- 
laxation and A^hb heatbath two-step updatings, 
changing ^oR and tuB links respectively. At a 
lattice spacing a ^ 0.17 fm and correlation length 
m-j^^rQ = 2.0 we simulated the HYP algorithm on 
a 8^ X 24 lattice with (^qr = 128, A^oR = 160) 
and (^hb = 200, A'hb — 80) with an acceptance 
of ~ 20%. If we run the standard thin link ac- 
tion with HMC at approximately matched physi- 
cal parameters, we obtain comparable autocorre- 
lations in the observables when the measurements 
are separated by a trajectory of unit time length. 
We can then compare the time costs of the two 
algorithms which are about 7 times larger for the 
HYP algorithm. 

With the new algorithm the quark mass can 
be lowered to values impractical with standard 
fermionic algorithms. A similar test as described 
above at m-^^ro = 1.6 shows that the HYP algo- 
rithm is about 3 times slower than the standard 
thin link HMC. 

A drawback of our algorithm is that as the lat- 
tice volume increases the numbers of the updated 
links ^oR and inB have to be kept unchanged, 
consequently the numbers of updatings Aqr and 
A^HB have to be scaled with the volume to keep 
the autocorrelation times unchanged. On the 
other hand as the continuum limit is approached 
the numbers of links toR and tuB which can be 
effectively updated scale: the physical volume of 
the updated region is constant. 

5. Uf ~ 2 HYP staggered fermions 

We consider eq. (|l|) with 
gtQ = (AftM)i/2, (12) 
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where M is the four flavor staggered fermion ma- 
trix. Perturbative arguments indicate that eq. 
(^2|) is equivalent in the continuum to a non-local 
theory of two flavors of quarks H . 



To evaluate the square root in eq. (|12|) we use 
a polynomial approximation. We write x^^^ = 
a;a;~i/2 = xP^i/2{x) with 

n 

P-i/2{x)= lim P^^IJx)= lim ycf)a;\(13) 
The coefficients Ci are found by minimizing |^ 

(14) 



I — I dx{x 
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For n even the roots of -P1"/2 come in complex 
conjugate pairs, we can write 



- „(») 



(»)t 



where 9["\/2) (■^) ~ 
allows us to identify 



Q = hm Mq 



(n) 

(-1/2) 



(M^M) . 



(15) 



). This 



(16) 



We note that the spectrum of M — 4(am)^ — 
Z?^ is bound from below exactly by 4(aTO)^ and 
the maximum eigenvalue fluctuates around the 
free field value 16 -I- 4(am)^. We set A = 1 in eq. 
( p^ and then rescale the polynomial to use it in 
the range x G (0, 18). 

In the accept/reject step eq. (^) the stochastic 
estimator exp AS* is evaluated as 



(18) 



where Af = M{V'), etc. and i? is a random 
Gaussian vector. The systematic error due to the 
finite order n of the polynomial is A^"^ = -Pi"/2 ~ 
P-i/2- It is possible to improve eq. (|l7|) such that 

(AS'("))i,„p = AS^"") + 0((A("))2) . (19) 

A similar improvement is possible for the vector 
^ in eq. ( p^ ) but we decided to take a higher 
order m = 128 of the polynomial. In the actual 
computation in order to increase the acceptance 
probability of the stochastic estimator we split 
the fermion matrix M = Mr A like in the Uf = 4 



case, eqs. ^ ^ 
In the evaluation of the polynomial eq. ( |l5| ) the 
ordering of the roots is important for the reduc- 
tion of numerical round-off errors. We emphasize 
that the smeared links make the polynomial ap- 
proximation work with lower order polynomials. 
In order to test the accuracy of the polynomial 
approximation we looked at the quantities 



^p{n,n') _ p(n) _ p{n'} 



(20) 



where P^cc = min{l, exp(AS'("))iinp} is the ac- 
ceptance probability using nth order polynomial, 
comparing n ~ 32, 64, 128 with n' ~ 256. We 
found that APacc = ©(lO^^) for n = 32 and de- 
creases to 0(10""*) for n = 128. Moreover APacc 
does not increase significantly as the mass de- 
creases or the lattice volume increases or the lat- 
tice spacing decreases. To answer the question of 
which accuracy is needed, we performed simula- 
tions with different values n = 32, 64, 128 but all 
other parameters fixed. First numerical results 
for the chiral condensate show agreement within 
the statistical errors for all orders n of the polyno- 
mial but much longer runs are needed to identify 
the systematic errors of a low order polynomial. 
The properties of the polynomial approximation 
have been studied extensively in and we plan 
to follow some of those methods in the future. 
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